AD-A081  18ft 

nwri  A9FIFTrn 


CHICAGO  UNI V  ILL  DEPT  OF  6E0PHYSICAL  SCIENCES  F/6  18/1 

reconstructive  tomography:  an  inverse  PROBLEM. (U) 

DEC  79  C  J  LEVINE  N00014-76-C-0039 

NL 


END 


DATE 

FUMED 

3  80 

DOC 


MICROCOPY  RESOLUTION  TEST  CHART 

NATIONAL  BUREAU  OF  STANDARDS  1963-A 


THE  UNIVERSITY  OF  CHICAGO 


RECONSTRUCTIVE  TOMOGRAPHY 


AN  INVERSE  PROBLEM 


A  DISSERTATION  SUBMITTED  TO 


THE  FACULTY  OF  THE  DIVISION  OF  THE  PHYSICAL  SCIENCES 


IN  CANDIDACY  FOR  THE  DEGREE  OF 


MASTER  OF  SCIENCE 


DEPARTMENT  OF  GEOPHYSICAL  SCIENCES 


CATHERINE  J.  LeVINE 


CHICAGO,  ILLINOIS 


DECfit®fi*r  *979 


- — n  ry 

‘  J  3 

; J "  1  'f 

Vs  <  21 

\f . 

U  to 

8a 

/ 

\ _ , 

TABLE  OF  CONTENTS 

Page 

ACKNOWLEDGEMENTS  .  111 

LIST  OF  ILLUSTRATIONS .  iv 

INTRODUCTION  .  1 

Chapter 

I.  EXISTENCE  AND  UNIQUENESS  OF  SOLUTIONS  .  9 

II.  RECONSTRUCTION  METHODS .  15 

III.  STABILITY  OF  SOLUTIONS . 32 

IV.  PROSPECTS  FOR  FUTURE  RESEARCH  .  35 


REFERENCES 


37 


ACKNOWLEDGEMENTS 


The  author  would  like  to  thank  Professor  Victor  Barcilon  for 
suggesting  this  paper,  and  for  encouragement  and  financial  support  during 
its  preparation. 

The  author  also  thanks  Dipl.  Geol.  Juliane  Fenner  for  translating 
the  paper  by  J.  Radon. 

The  author  also  thanks  Marilyn  Bowie  for  typing  the  manuscript. 

The  author  would  like  to  express  special  appreciation  to  Stephen 
Daly,  David  Thorn  and  Robert  Haskins  for  their  continued  encouragement 
ar.d  moral  support. 

This  research  was ^sppported  by  the  Office  of  Naval  Research  under 
contract  N00014-76-C-0034 .  Reproduction  in  whole  or  in  part  is  permitted 
for  any  purpose  of  the  United  States  Government. 


iii 


LIST  OF  ILLUSTRATIONS 


Figure  Page 

1  Geometry  of  line  integral  data  for  parallel 

mode  tonography .  41 

2  Geometry  of  fan  bean  scanning .  42 

3  Star-shaped  artifact  in  back  projection  method .  43 

4  Convolution  or  filtered  back  projection  method .  44 

5  Spatial  discretization  of  reconstruction  domain 

and  line  integral  data . . .  45 


iv 


INTRODUCTION 


■i 

The  process  of  recovering  the  three-dimensional  structure  of  an 

object  by  reconstructing  successive  cross  sections  orthogonal  to  a 

common  axis  is  known  as  reconstructive  tomography.  Each  cross  section 

of  the  object  is  essentially  a  function  f  defined  on  a  two-dimensional 

domain  C.  This  unknown  function  may  be  specified  using  as  data  integrals 

of  f  along  straight  lines  traversing  C.  Data  of  this  sort  represent 

cumulative  measures  of  f  along  certain  paths  rather  than  actual  values. 

of  f  at  specific  points .  .  Suppose  f  vanishes  outside  the  unit  disk  D: 

?  2  ~ 

x"  -r  y  $  1;  since  any  finite  domain  may  be  enclosed  by  a  disk  this  is 
not  a  restriction  on  C.  Let  L(t,9):  x  cos  e  +  y  sin  9  =  t  denote  the 
straight  line  which  makes  an  angle  9  with  the  positive  y  axis  and  whose 
perpendicular  distance  to  the  origin  is  t,  as  shown  in  Fig.  1.  The  integral 
of  f  along  L(t,  9)  is 

p(t,  C)  =  /f(?)  ds  (1) 

L(t , 9) 

A  A 

where  r  =  (x,  y)  and  £  is  the  unit  vector  £  =  (cos9 ,  sin0) .  The  original 
coordinates  x,y  are  related  to  the  rotated  coordinates  t,s  by 

x  =  t  cose  -  s  sin9  ^2) 

y  =  t  sin9  +  s  cos9 . 


Equivalently,  the  integral  of  f  along  L(t,9)  may  be  expressed  as 


where  r*£  ■  x  cos8  +  y  sin9  and  <5(r)  is  the  Dirac  delta  function.  The 
set  of  line  integrals  for  various  values  of  t,  |t|  s  1,  at  a  given  angle 
6  is  known  as  a  scan  of  f.  When,  for  a  given  9,  the  integral  values  are 
specified  for  all  t,  |t|  j  1,  the  scan  will  be  called  a  projection  of  f 
onto  the  line  orthogonal  to  the  direction  9.  In  much  of  the  following 
discussion  the  line  integral  data  will  be  assumed  to  be  projections  of  f, 
that  is,  continuous  in  t.  A  complete  data  set  includes  projections  for 
all  angles  9e[0,2T:).  This  problem  has  applications  in  a  number  of  fields. 

The  reconstruction  of  the  human  body  has  received  much  attention 
in  the  last  decade.  Computerized  Axial  Tomography  (CAT)  transmission 
scanners  use  measurements  of  transmitted  radiation  such  as  x-rays,  gamma 
rays  or  protons  to  determine  the  corresponding  attenuation  coefficient 
and,  by  inference,  the  density  of  a  cross  section  (Cormack  1963,  1964; 

Che  st  at.  1974).  The  data  represent  the  total  attenuation  of  a  beam  of 
radiation  due  to  its  passage  through  an  object.  Emission  scanners  recon¬ 
struct  the  isotope  distribution  of  a  cross  section  from  measurements  of 
radionuclide  emissions  (Budinger  and  Gullberg  1.974b,  Kuhl  and  Edwards 
1968,  Kuhl  et  at.  1973).  Most  scanners  operate  in  the  parallel  mode  in 
which  the  line  integral  data  are  obtained  by  taking  measurements  at  succes¬ 
sive  positions  of  the  detector  as  it  is  translated  parallel  to  the  t  axis 
as  shown  in  Fig.  1.  In  the  diverging  beam  or  fan  beam  scanners  the  object 
is  scanned  by  radiation  diverging  from  a  point  source  S  as  shown  in  Fig.  2. 
With  either  configuration  the  source-detector  apparatus  is  rotated  about 
the  object  or,  equivalently,  the  object  is  rotated  within  the  apparatus 
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to  produce  a  set  of  line  integral  scans.  Since  the  density  or  the  isotope 
concentration  of  abnormal  tissue  differs  from  that  of  healthy  tissue,  such 
reconstructions  are  a  useful  diagnostic  tool.  They  are  especially  important 
in  the  detection  of  brain  tumors,  a  task  for  which  conventional  x-rays  are 
of  little  use. 

The  three-dimensional  structure  of  bionolecules  may  be  determined 
from  electron  micrographs.  Rotation  of  the  specimen  within  the  micro¬ 
scope  results  in  a  set  of  two-dimensional  maps  of  the  total  attenuation 
of  the  electron  beam  as  it  passes  through  the  particle  in  different  orien¬ 
tations.  The  three-dimensional  optical  density  of  the  specimen  may  be  re¬ 
constructed  directly  from  this  set  of  two-dimensional  projections  (Hart 
1968;  DeRosier  and  Klug  1968;  Crowther  et  at.  1970) .  In  an  alternative 
method  (Gordon  et  at.  1970;  Vainshtein  1970,  1971),  each  micrograph  is 
regarded  as  a  set  of  parallel  lines;  each  line  is  a  measure  of  the  attenu¬ 
ation  of  the  beam  due  to  passage  through  a  particular  cross  section  of 
the  particle.  Each  cross  section  is  reconstructed  from  the  set  of  these 
one-dimensional  projections,  one  from  each  micrograph,  representing  dif¬ 
ferent  orientations  of  the  particle.  The  three-dimensional  density  struc¬ 
ture  is  then  recovered  by  stacking  the  cross  sections. 

In  radio  astronomy  the  two-dimensional  brightness  distribution  of 
radio  sources  or  the  Sun  or  Moon  are  inferred  from  measurements  of  the 
intensity  of  radiation  received  by  radio  telescopes  (Bracewell  and  Roberts 
1954,  Bracewell  and  Riddle  1967).  The  telescope  collects  radiation  along 
a  narrow  strip  which  may  be  idealized  as  a  line.  The  amount  of  radiation 
received  for  a  given  orientation  of  the  antenna  is  a  relative  measure  of 
the  total  radiation  emitted  along  a  similarly-positioned  line  across  the 
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surface  of  the  source.  Adjustment  of  the  telescope  or  the  rotation  of 
the  Earth  results  in  a  scan  of  lines  at  the  same  angle.  A  set  of  scans 
is  obtained  by  changing  the  orientation  of  the  antenna.  Records  of  a 
lunar  occultation  of  a  star  will  produce  similar  scans  and  the  star's 
brightness  distribution  may  be  reconstructed  from  data  for  occultations 
occurring  from  different  directions. 

The  use  of  travel  time  data  for  waves  propagating  through  a  medium 
allows  properties  of  the  medium  to  be  inferred.  Such  data  represent  the 
integrals  of  the  slowness,  the  reciprocal  of  the  velocity,  along  the 
paths  of  the  waves.  In  seismology  travel  times  traditionally  have  been 
used  to  determine  the  vertical  seismic  velocity,  and  hence  density,  struc¬ 
ture  of  the  earth.  Recently  Aki  et  at.  (1977)  developed  a  method  for 
finding  the  three-dimensional  slowness  structure.  They  assumed  an  ini¬ 
tially  layered  earth  model  and  divided  each  layer  into  blocks  within 
which  a  perturbation  to  the  initial  slowness  model  is  assumed  to  be  con¬ 
stant.  The  travel  tine  residual,  the  difference  between  the  observed 
travel  tine  and  that  calculated  for  the  initial  model,  for  each  seismic 
ray  received  is  then  the  sum  of  the  slowness  parameters  for  each  block 
through  which  the  ray  has  passed.  With  this  method  the  problem  is  solved 
in  three  dimensions.  Munk  and  Wunsch  (1979)  treated  the  problem  of 
determining  mesoscale  fluctuations  in  the  density  of  the  ocean  from 
travel  times  for  acoustic  waves  in  a  similar  manner.  They  used  the 
block  structure  of  Aki  et  at.  to  reconstruct  the  sound  speed  perturbation 
within  either  horizontal  or  vertical  cross  sections  and  stacked  the  cross 
sections  to  obtain  the  three  dimensional  structure.  Their  data,  based  on 
an  array  of  several  acoustic  sources  and  receivers,  was  composed  of  all 
the  resolvable  waves  arriving  at  each  receiver  along  all  possible  ray  paths. 


Some  of  these  applications  involve  the  reconstruction  of  a  three- 
dimensional  function  from  a  set  of  two-dimensional  projections  or  the  re¬ 
construction  of  two-dimensional  functions  from  line  integral  data  rather 
than  the  reconstruction  of  a  three-dimensional  function  by  stacking  a 
sequence  of  cross  sections  obtained  from  line- integral  data.  These  appli¬ 
cations  are  not  considered  tomography  in  the  strict  sense  but  techniques 
devised  for  the  solution  of  these  problems  have,  been  applied  to  the  tomo¬ 
graphy  problem  as  well. 

Each  of  the  problems  described  above  may  be  characterized  by  the 
formulation  in  Eq.  (3).  This  equation,  restated  in  canonical  form,  is 

p(t,C)  =  If  =  f  G(r,5;t)  f (r)  d2r  (4) 

A  A 

where  G(r,£;t)  =  5(r*£  -  t)  and  t  is  a  parameter.  Since  the  operator  L 
is  linear  the  problem  of  solving  this  equation  for  f  may  be  approached 
using  well-developed  techniques  in  inverse  theory.  The  prediction  of  the 

^  a 

data  p(t,£)  when  f(r)  is  known  is  the  direct  or  forward  problem.  In  general 
the  direct  problem  involves  the  formulation  of  a  model  relating  some  observ¬ 
able  property  of  a  function  to  the  function  itself.  The  appropriateness 
of  the  model  is  evaluated  by  comparison  of  calculated  data  with  observed 
data.  The  tomography  problem  involving  transmitted  radiation,  for  example, 
is  based  on  the  observed  exponential  attenuation  of  radiation  due  to 
passage  through  matter.  The  intensity  I(t,£)  of  a  transmitted  beam  of 
radiation  is  known  to  be  related  to  the  linear  attenuation  coefficient 
u(r)  of  the  medium  through  which  it  has  passed  by 


in  which  I  is  the  intensity  of  the  incident  ratiation.  The  data 

p(t,C)  =*  -  ln[I(t,0/I0]  (6) 

are  thus  related  to  u(r)  by  Eq.  (1).  A  similar  formulation  of  a  model 
which  relates  a  function  f  to  a  measurable  property  of  f  is  more  difficult 
because  the  unknown  is  a  function  rather  than  a  collection  of  numbers. 

For  this  reason  solution  of  the  inverse  problem  is  generally  attempted 
by  inverting  the  operator  L,  resulting  in 

f  (r)=Ir-1[p(t,C)  ]  (7) 

where  L  ^  is  the  operator  inverse  to  L,  In  most  problems  this  is  not 
an  easy  task  either.  In  addition  to  the  construction  of  a  solution  which 
satisfies  Eq.  (4)  for  a  specified  data  set,  inverse  theory  addresses  such 
questions  as  the  existence,  uniqueness  and  stability  of  possible  solutions. 

The  purpose  of  this  paper  is  to  review  the  literature  on  recon¬ 
structive  tomography  within  the  framework  of  inverse  theory  to  be  described 
below.  The  concepts  of  existence,  uniqueness,  construction  and  stability 
of  solutions  to  Eq.  (4)  will  be  explored  and  discussed  in  terms  of  their 
application  to  tomography. 

Proof  of  the  existence  of  a  solution  to  an  inverse  problem  ensures 
that  there  is  a  solution  which  is  consistent  with  the  problem  as  posed;  it 
does  not  necessarily  imply  that  a  solution  has  been  found.  Often  the 
assumptions  incorporated  into  the  formulation  of  a  problem  impose  con¬ 
straints  on  the  data.  If  these  consistency  conditions  are  violated  there 
is  no  solution  which  satisfies  both  the  data  and  the  original  assumptions. 

In  the  tomography  problem,  for  example,  if  the  solution  is  assumed  to  be 
radially  symmetric  its  scans  should  be  the  same  for  all  angles.  If  the 
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measured  scans  are  not  in  fact  independent  of  9,  the  assumption  of  radial 
symmetry  is  invalid.  The  consistency  of  a  data  set  with  such  constraints 
does  not,  however,  guarantee  the  existence  of  a  solution.  In  many  prob¬ 
lems  a  solution  which  is  consistent  with  the  data  is  assumed  to  exist. 

If  a  solution  to  an  inverse  problem  does  exist,  it  may  not  be 
the  only  solution  which  satisfies  the  given  data.  Possible  sources  of 
such  non-uniqueness  include  contamination  of  the  data  by  errors  or,  in 
the  case  of  perfect  data,  the  finiteness  of  the  data  set.  The  non¬ 
uniqueness  associated  with  a  partial  data  set  may  be  motivated  in  terms 
of  the  tomography  problem  as  follows:  There  are  an  infinite  number  of 
functions  which  satisfy  a  data  set  composed  of  only  one  scan.  A  given 
line  integral  value  p(t,  £)  might  be  produced,  for  instance,  by  a  point 
source  at  any  of  the  points  along  L(t,  9);  there  is  no  way  to  distinguish 
between  these  possible  solutions.  If  a  second  scan  is  introduced  the 
number  of  solutions  which  produce  both  scans,  while  still  infinite,  is 
reduced.  Only  when  scans  are  given  for  all  angles  9e(0,27r)  is  a  unique 
solution  possible.  For  some  inverse  problems  even  a  complete,  error- 
free  data  set  is  insufficient  to  resolve  the  non-uniqueness.  See  Backus 
(1S70)  for  an  example  of  such  inherent  non-uniqueness.  Even  though  a 
solution  is  completely  undetermined  by  a  finite  data  set,  some  information 
may  be  extracted.  Since  obtaining  a  unique  solution  is  generally  not 
possible  a  reasonable  alternative  is  to  attempt  to  characterize  the  class 
of  solutions  which  satisfy  the  given  data.  Bounds  on  some  property,  such 
as  the  minimum  value,  of  the  solution  may  be  established  by  selecting 
from  this  class  that  solution  which  extremizes  the  property  of  interest. 
Such  extremal  solutions  are  unique  (Parker  1972,  1975;  Sabatier  1977a, 


1977b)  . 
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An  inverse  problem  is  considered  well-posed  or  stable  if  for  each 
set  of  data  a  solution  exists,  is  unique  and  depends  continuously  on  the 
data.  Continuous  dependence  on  the  data  implies  loosely  that  a  small 
change  in  the  data  will  produce  a  small  change  in  the  solution.  This 
means  that  an  error  in  the  data  will  not  cause  the  solution  to  deviate 
substantially  from  the  solution  which  would  be  obtained  if  the  data  were 
perfect.  Linear  inverse  problems  on  a  finite  domain  described  by  an 
integral  equation  of  the  form  (4)  are  generally  unstable  if  the  kernel 
G  is  continuous. 

There  is  often  more  than  one  method  of  constructing  a  solution 
to  an  inverse  problem.  This  is  especially  true  for  the  tomography  prob¬ 
lem.  Construction  procedures  are  generally  judged  on  their  accuracy  and 
the  computation  time  they  require. 

This  paper  will  focus  on  reconstructive  tomography  as  an  inverse 
problem.  There  are  other  review’s  which  emphasize  different  aspects  of 
this  subject.  3rooks  and  DiChiro  (1976)  reviewed  several  reconstruction 
methods  and  discussed  experimental  and  diagnostic  limitations  of  tonography. 
Gordon  and  Herman  (1974)  presented  a  comparative  survey  of  different  re¬ 
construction  techniques.  Reviews  by  Smith  et  al.  (1977)  and  Shepp  and 
Kruskal  (1978)  described  contributions  to  the  mathematical  theory  of 
tomography. 


CHAPTER  I 


EXISTENCE  AND  UNIQUENESS  OF  SOLUTIONS 

The  operator  L  defined  in  Eq.  (4)  is  known  as  the  Radon  transform 

after  Johann  Radon  who  originally  solved  the  reconstruct  ion  problem 

(Radon  1917).  Radon  formed  an  average  p(r;q)  of  the  line  integral 

-2 

values  at  a  distance  q  from  r  =  (x,y) ,  weighted  it  by  a  ,  and  summed 
over  all  possible  values  of  q  to  obtain  the  value  of  the  source  function 
at  r.  His  inversion  formula  is 


1  lim  P(^;e) 
*  e-*0  e 


(8) 


(9) 


where 


P(r;q) 


2tt 

p(r*C  -I-  q,C) 


d9. 


(10) 


_2 

Weighting  the  average  by  q  emphasizes  the  contribution  of  line  integrals 
near  r  and  lessens  that  of  those  at  a  distance.  Radon  assumed  that  the 
line  integral  scans  p(t,E)  are  continuous  in  t,  are  specified  for  all 
angles  9e[0,2r),  have  continuous  derivatives  in  t  and  8,  and  that  both 
the  scans  and  their  derivatives  satisfy  certain  conditions  on  their  asymp¬ 
totic  behavior.  Under  these  assumptions  he  proved  that  the  function  £(r) 
obtained  by  this  inversion  formula  exists,  is  real  and  continuous,  and  is 
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unique.  Thus  in  the  tomography  problem  the  non-uniqueness  associated  with 
a  partial  data  set  is  resolved  when  a  complete  data  set  is  given. 

Radon  also  generalized  the  inversion  to  n  dimensions  using  as 
data  integrals  over  (n-1) -dimensional  hyperplanes.  If  n  is  even  the  re¬ 
construction  of  f  at  a  point  r  involves  all  the  line  integral  data  whereas 
if  n  is  odd  only  those  integrals  along  lines  L(t,9)  in  a  small  neighborhood 
of  r  are  required  (Gel1 f and  et  at.  1966).  The  reconstructions  discussed  in 
this  paper  will  be  two-  and  three-dimensional. 

In  order  for  a  set  of  functions  {g(t,£):  9ef0,2ir)}  to  represent 

scans  of  a  solution  f  the  functions  g(t,£)  must  satisfy  certain  consistency 

conditions  imposed  by  the  form  of  the  operator  L.  One  such  condition  is 

that  g(t,5)  =  g(-t,-')>  that  is,  that  the  integral  of  f  along  any  line 

L(t,9)  must  be  independent  of  the  direction  of  integration.  Secondly, 

the  functions  g(t,£)  should  belong  to  the  same  space,  with  a  different 

domain  of  definition,  as  the  function  f  to  be  reconstructed.  For  example, 

a  function  g(t,£)  representing  a  scan  of  a  function  which  vanishes  outside 

the  unit  disk  must  vanish  for  jt|  >1.  A  final  condition  is  that 
r*  *  k 

I  g(t,Ot  dt  be  a  polynomial  in  £  of  degree  less  than  or  equal  to  any 
positive  integer  k.  This  condition  is  a  consequence  of  the  equivalence 


£ 


(LfXt.O  t*  dt  = 


f  f  (r)  (r-Ok  d2r 
*/D 


(ID 


in  which  the  right  side  is  a  polynomial  in  5  of  degree  less  than  or  equal 
to  k.  These  conditions  are  both  necessary  and  sufficient  for  the  exis¬ 
tence  of  a  function  f  such  that  g  =  Lf  (Ludwig  1966).  Other  consistency 
conditions  arise  in  connection  with  certain  reconstruction  methods  and 


will  be  mentioned  as  these  methods  are  described. 
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The  uniqueness  result  obtaiued  by  Radon  is  in  general  valid  only 
when  scans  at  all  angles  8e[0»2ir)  are  available.  One  exception  occurs 
when  the  unknown  function  is  known  to  be  radially  symmetric.  In  this 
case  one  scan  which  is  continuous  in  t  is  equivalent  to  the  complete  set 
of  scans  and  f(r)  is  determined  uniquely  (Bracewell  1956).  A  second  ex¬ 
ception  provides  a  measure  of  the  resolution  which  may  be  achieved  in  a 
reconstruction  based  on  a  finite  number  of  scans.  Suppose  a  function 
f(r)  is  constant  on  square  grid  elements  with  side  length  b  and  vanishes 
for  Jx|  >  a,  |y|  >  a,  where  a  is  a  multiple  of  b.  Then  f(r)  may  be  re¬ 
constructed  uniquely  from  a  set  of  scans  at  angles  0  ^ ,  j-1,  ....  2a/b, 
such  chat  0  <  <  tan  ^  1/ ( j  — 1)  <  u/2  (Frieder  and  Herman  1971).  This 

restriction  on  the  spacing  of  the  angles,  which  rules  out  a  unique  recon¬ 
struction  if  the  angles  are  equally  spaced,  is  actually  an  advantage  in 
electron  microscopy  where  the  range  of  available  angles  is  generally  very 
narrow.  The  scans  may  be  sampled  discretely  in  t  if  the  sampling  interval 
is  lass  than  /~2  b.  Any  feature  smaller  than  b  which  appears  in  a  recon¬ 
struction  must  be  an  artifact  of  the  reconstruction  procedure. 

When  line  integral  scans  at  a  finite  number  of  angles  9^,j=l,  ...,  N, 
are  given  the  reconstruction  problem  may  be  stated  as 

pCt.Cj)  -  £jf  =  f  Hr)  6(r-ij  -  t)  d2r  j  -  1,  ....  N  (12) 

•'d 


where  5  =  (cos9^,  sinB^).  In  the  absence  of  a  priori  knowledge  concerning 

N 

the  symmetry  properties  of  the  unknown  f,  the  data  set  {p(t,£j)J  is 

not  sufficient  to  specify  f  uniquely  (Solomon  1976).  There  will  always 


be  functions  which  integrate  to  zero  in  the  directions  0 ^ .  Solomon  proved 
that  any  function  in  the  null  space  of  L is  orthogonal  to  all  functions 
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which  are  constant  along  the  lines  L(t,6^),  |t|  S  1,  and  gave  a  formula 
for  the  angles  between  the  null  spaces  in  terms  of  the  angles  0^ .  Bracewell 
and  Roberts  (1954)  also  showed  that  there  are  "invisible  distributions" 
which  cannot  be  resolved  by  a  finite  number  of  projections.  They  defined 
a  principal  solution  to  be  that  solution  which  satisfies  the  given  data 
and  for  which  the  null  spaces  are  assumed  to  be  zero. 

Use  of  a  sampling  theorem  for  Fourier  transforms  allows  this 
principal  solution  to  be  determined  uniquely  from  a  finite  number  of 
scans.  A  function  may  be  completely  recovered  by  a  discrete  sampling 
of  its  two-dimensional  Fourier  transform  if  the  transform  vanishes  out¬ 
side  a  circle  of  radius  Rq,  known  as  the  cut-off  frequency.  If  the  source 

function  f(r)  is  confined  within  a  circle  of  radius  r  ,  the  sampling  inter- 
~  o 

val  should  be  less  than  rQ  ^  (Bracewell  and  Riddle  1967).  Knowledge  of 
the  scans  of  f  in  the  directions  9  ^ ,  j=l,  . ..,  N,  is  equivalent  to  knowl¬ 
edge  of  the  two-dimensional  Fourier  transform  of  f  along  radial  lines  at 
the  same  angles,  as  will  be  detailed  below.  If  the  number  of  angles  is 

chosen  so  that  the  spacing  of  the  radial  lines  near  R=R  is  within  the 

o 

sampling  interval  rQ  then  a  principal  solution  for  f  is  completely 
determined,  up  to  a  component  corresponding  to  spatial  frequencies  greater 
than  Rq.  This  idea,  developed  in  radio  astronomy,  is  based  on  the  prop¬ 
erties  of  radio  telescopes  which  cannot  resolve  frequencies  beyond  some 
cut-off  frequency  Rq.  Although  the  Fourier  transform  of  an  actual  distri¬ 
bution  may  not  vanish  for  R  >  Rq,  there  is  no  loss  of  information  involved 
in  making  this  assumption  since  such  frequencies  cannot  he  measured.  The 
principal  solution  represents  all  the  information  which  may  be  obtained 


in  this  manner  from  a  finite  set  of  scans.  If  the  amplitude  of  the 
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transform  is  small  near  the  cut-off  frequency,  implying  a  nearly  continuous 
transition  to  zero,  the  principal  solution  will  be  a  good  approximation  to 
the  true  solution  (Bracewell  and  Roberts  1954) .  Since  one  of  the  necessary 
conditions  for  the  existence  of  a  solution  to  the  tomography  problem  is 

A 

that  p(t,£)  *  0  for  |t|  >1,  the  two-dimens ioiial  Fourier  transform  of  f, 
determined  from  the  one-dimensional  transform  of  p(t,£),  will,  in  most 
instances,  be  small  near  R=*l  and  for  R  >  1  (Bracewell  1965,  Ch.  18).  As 
a  result  the  principal  solution  obtained  from  a  finite  set  of  scans  will 
be  a  reasonable  approximation  to  the  actual  solution.  If  the  amplitude 
of  the  transform  is  not  small  near  R=1  the  principal  solution  will  not 
be  a  good  approximation  to  the  true  solution;  the  discontinuity  in  the 
transform  produces  oscillations  in  the  solution. 

Logan  (1975)  expressed  this  result  more  quantitatively.  Let 
e(r)  =  f(r)  -  f.T(r)  where  f(r)  is  the  true  solution  and  f,.(r)  an  approxi- 
nation  based  on  a  set  of  N  scans.  Let  E(R)  and  F(R)  denote  the  Fourier 
transforms  of  s(r)  and  f(r),  respectively.  If  F(R)  is  concentrated  within 
a  circle  of  radius  R  ,  define  A„(R  )  as  the  fraction  of  the  total  energy 
of  E(R)  which  is  contained  within  the  same  circle.  This  cut-off  frequency 
Rc  is  the  same  as  that  mentioned  above  if  F(R)  actually  vanishes  for 
R  >  Ro;  otherwise  Rc  >  Rq.  If  little  of  the  energy  of  E(R)  is  contained 
within  Rc,  that  is,  if  ^N(RC)  is  small,  then  e  is  small  and  f  is  approxi¬ 
mately  reconstructed  by  N  scans.  If  A.,(R  )  is  near  one  then  e  has  a  large 

N  c 

component  on  the  null  space  of  the  operator  L  and  fN  is  a  poor  approxi¬ 
mation  to  f.  Logan  showed  that 

A  (N  +  ctNP)  =1  a  >  0,  p  >  1/2 
N 

(13) 

A^(N  -  aNP)  =0  a  >  0,  p  >  1/3. 


lim 

N-*=° 

lim 

N-"» 
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Thus  a  set  of  N  scans  can  resolve  a  function  f (r)  whose  Fourier  transform 

is  contained  within  a  circle  of  radius  R  which  is  somewhat  less  than  N. 

c 

Conversely,  if  the  Fourier  transform  of  f  is  known  to  be  contained  within 
a  circle  of  radius  Rc>  then  N,  somewhat  greater  than  R^,  projections  are 
needed  to  obtain  a  good  approximation  to  f. 

A  finite  data  set,  although  generally  insufficient  to  specify 
a  unique  solution  to  an  inverse  problem,  does  limit  the  class  of  possible 
solutions.  One  means  of  characterizing  this  reduced  class  of  viable 
solutions  is  through  extremal  solutions.  Parker  (1972)  presented  a 
variational  approach  for  finding  bounds  on  various  properties  of  possible 
model  solutions,  such  as  an  integral  of  the  model,  its  norm  or  its  maximum 
value.  In  his  terminology  (Parker  1975),  the  ideal  solution  is  that  solu¬ 
tion  which  extremizes  some  given  property  of  the  model.  For  the  recon¬ 
struction  problem  Logan  and  Shepp  (1975)  proved  the  uniqueness  of  the 
ideal  solution  having  the  minimum  norm  and  showed  it  to  be  of  the  form 


(r)  =  l  p  (r •£  ). 

j=l  3  3 


The  Pj's>  known  as  ridge  functions,  are  constant  along  the  lines  L(t,0^) 
and  vary  only  in  t.  They  gave  an  explicit  form  for  f(r)  when  the  angles 
0^  are  equally  spaced.  Marr  (1974)  developed  a  least  squares  polynomial 
solution  to  the  diverging  beam  tomography  problem.  He  reconstructed  a 
function  on  the  unit  disk  using  as  data  line  integrals  along  the  M(M-l)/2 
chords  which  connect  each  pair  of  M  equally  spaced  points  on  the  unit 
circle.  He  proved  that  within  the  space  of  polynomials  the  solution 
which  best  fits  the  data  in  a  least  squares  sense  is  a  unique  polynomial 
,  in  x  and  y  of  to.tal  degree  d  $  M-2. 


CHAPTER  II 


RECONSTRUCTION  METHODS 

Radon's  inversion  formula  is  impractical  for  implementation  because 
it  is  difficult  to  choose  a  discretization  so  that  there  are  enough  inte¬ 
gral  values  for  lines  at  a  distance  q  from  r  *  (x,y)  to  give  a  good 
average  (Shepp  and  Kruskal  1978).  As  a  result  many  methods  for  the 
approximate  reconstruction  of  functions  from  their  line  integrals  have 
been  devised.  In  most  methods  the  line  integral  data,  given  for  a  finite 
number  of  angles  S  ,  j=l,  ...,  N,  will  be  assumed  to  be  continuous  in  t. 

Summation  methods 

The  summation  methods  involve,  as  the  name  implies,  a  super¬ 
position  of  the  line  integral  data.  Hart  (1968)  used  a  photographic 
superposition  in  electron  microscopy.  The  general  shape  of  biological 
particles  is  inferred  from  the  bright  region  where  the  micrographs  over¬ 
lap.  A  second  type  of  summation  is  known  as  back  projection  (Vainshtein 
1970).  In  this  method,  also  originally  developed  in  electron  microscopy, 
the  line  integral  scans  are  regarded  as  the  projections  of  the  object  to 
be  reconstructed.  The  projections  at  different  angles  are  equidistant 
from  the  center  of  the  object  and  may  be  thought  of  as  tangents  to  a 
circle  K  concentric  to  and  enclosing  the  domain  D.  The  point  of  tangency 
of  each  projection  is  the  point  t=Q.  For  each  line  L(t,0^)  in  the  scan 

A 

at  the  angle  8^  the  corresponding  value  of  p(t,£^)  is  assigned  to  each 
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point  along  that  line  within  the  circle  K.  This  process  is  called  back 
projecting.  The  values  of  pCt.C^)  are  assigned  in  a  similar  manner;  at 
those  points  along  L(t,02)  to  which  a  value  was  assigned  during  the  back 
projection  of  p(t,£^)  the  appropriate  value  of  pCt.C^)  is  added  to  the 
previously  assigned  value.  This  procedure  is  continued  until  each  scan 
has  been  back  projected.  The  value  of  the  reconstructed  function  at  a 
point  r  is  then  the  sum  of  all  the  line  integral  values  p(t,£j)  corre¬ 
sponding  to  t's  for  which  L(t,9^)  intersects  r.  Kuhl  and  Edwards  (1968) 
adapted  this  method  slightly  by  averaging  the  line  integral  values  assigned 
to  each  point.  A  major  drawback  of  both  the  photographic  superposition 
and  Che  back  projection  methods  is  that  the  reconstructed  function  does 
not  satisfy  the  original  line  integral  data.  Kuhl  et  al.  (1973)  solved 
this  problem,  in  an  approach  called  orthogonal  tangent  correction,  by 
back  projecting  only  two  orthogonal  scans  at  a  time  and  correcting  the 
solution  at  each  step  so  it  satisfies  all  the  scans  used  in  previous  steps. 

A  second  disadvantage  of  the  back  projection  method  is  that  it  produces  a 
star-shaped  pattern,  shown  in  Fig.  3,  which  arises  when  the  back  projections 
of  two  or  more  scans  overlap  in  regions  outside  the  domain  D.  This  over¬ 
lapping  outside  D  makes  it  difficult  to  ascertain  the  actual  extent  of 
the  reconstructed  solution.  An  even  more  important  limitation  of  this 
method  is  that  the  function  it  reconstructs  is  actually  a  blurred  version 
of  the  true  solution.  Remedies  for  the  last  two  problems  will  be  discussed 
below. 

Fourier  transform  method 

Due  to  the  nature  of  the  operator  L,  the  one-dimensional  Fourier 
transform  of  a  line  integral  scan  p(t,£)  of  f  is  equivalent  to  the  corre- 
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sponding  central  section  of  the  two-dimensional  Fourier  transform  of  f. 
The  central  section  referred  to  is  the  line  passing  through  the  origin 
and  making  an  angle  6  with  the  horizontal  axis  in  the  transform  plane. 
Let  the  one-dimensional  transform  of  p(t,£)  be  defined  as 


P(T 


,,  ,f 

•'•a 


p(t,£)  exp  C-2ititT)  dt 


(15) 


I 


It/ 6  (r  > 
Jn 


=  I  exp  (-2iritT)  dt|6  (r*^  -  t)  f(r)  d  r 
'D 


(16) 


Interchanging  the  order  of  integration  and  performing  the  integration 
eve-  t  results  in 


P(T,£) 


exp  (2ir  ir'CDd^r  . 


(17) 


If  the  two-dimensional  transform  of  f  is  defined  as 


F(R) 


exp(2irir-R)d2r 


(18) 


Eq.  (17)  may  be  restated  as 

P(T.O  “  F(t£).  (19) 

In  the  Fourier  transform  method  F(T£)  is  constructed  by  transforming  the 
line  integral  scans  and  the  solution  is  obtained  by  inverting  F(TC) 

(Crowther  et  al.  1970).  If  scans  at  all  angles  0e[O,2ir)  are  available 
the  Fourier  transform  of  f  is  completely  specified  and  the  solution  f 
is  unique.  Even  when  the  set  of  scans  is  finite  this  method  can  provide 
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a  good  approximation  to  f.  As  discussed  in  the  last  section  the  principal 

solution  of  Bracewell  and  Roberts  (1954)  may  be  determined  uniquely  by  a 

discrete  sampling  of  its  Fourier  transform.  If  the  number  of  scans  is 

large  enough  to  permit  an  adequate  sampling  of  the  transforms  P(T,£) 

this  principal  solution,  generally  a  good  approximation  to  the  true 

solution,  is  recovered  uniquely.  Bracewell  and  Riddle  (1967)  showed  that 

in  order  to  reconstruct  an  object  of  diameter  D  there  should  be  N  **  irD  R 

equally  spaced  scans.  Other  authors  (Gilbert  1972b,  Logan  1975)  have  also 

discussed  this  trade-off  between  the  resolution  d  =  R  ^  and  the  number  of 

o 

scans.  DeRosier  and  Klug  (1968)  and  Vainshtein  (1970)  noted  that  symmetry 
considerations  reduce  the  number  of  scans  needed  to  obtain  a  certain  reso¬ 
lution.  Bracewell  and  Roberts  (1954)  suggested  possible  modifications  to 
be  used  when  the  amplitude  is  not  small  near  the  cut-off  frequency  although 
such  corrections  reduce  the  resolution  of  the  solution.  Although  the 
Fourier  transform  method  produces  good  reconstructions,  two-dimensional 
inversions  are  computationally  very  time  consuming;  this  is  the  primary 
drawback  of  this  method.  If  F(R)  is  expressed  in  polar  coordinates  a  one- 
dimer.sional  Fourier-Bessel  inversion  may  be  used  (Crowther  et  at.  1970); 
this  inversion  requires  less  computing  time  but  reduces  the  resolution. 

Gilbert  (1972b)  developed  a  reconstruction  method  based  on  the 
equivalence  (19)  which  does  not  require  a  two-dimensional  inversion.  This 
method  makes  use  of  the  convolution  theorem  for  Fourier  transforms  which 
states  that  the  Fourier  transform  of  the  convolution  of  two  functions  is 
the  product  of  their  transforms  (Bracewell  1965,  Ch.  6).  If  the  transform 
F(T£)  used  in  a  numerical  inversion  is  regarded  as  the  product  of  the  con¬ 
tinuous  transform  and  a  sampling  function  S(T,9), 
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which  samples  the  transform  with  unit  weight  at  the  origin  and  along 
equally  spaced  radial  lines  corresponding  to  the  angles  of  the  scans. 

Using  this  sampling  function  Gilbert  derived  the  following  relation 
between  the  function  reconstructed  using  the  back  projection  method 
(Vainshtein  1970)  and  the  true  solution: 

f.  (r)  =  f (r)  *  s(t,9) .  (23) 

bp  ~ 


In  the  limit  as  N-**> 

VS>  ■  f<5>  *  IVT  ■  <“> 

a  result  also  obtained  by  Vainshtein  (1971).  This  formulation  illus¬ 
trates  the  primary  limitation  of  the  back  projection  method,  that  the 
reconstructed  function  is  in  fact  the  actual  solution  convolved  with 
t  In  terms  of  Gilbert's  method  this  is  a  result  of  sampling  the 
transform  plane  with  equal  weight  at  all  radii  without  compensating 
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for  the  overlapping  of  the  scans  near  t=0.  The  sampling  function 
S’(T,0)  =  T  S(T,6)  produces  better  results.  Since  the  error  introduced 
by  the  convolution  of  the  solution  with  S’(T,0)  is  of  the  same  order 
as  the  truncation  error  in  a  two-dimensional  Fourier  inversion  the  two 
approaches  are  equivalent. 

Weighting  by  the  transform  variable  T  was  also  proposed  by 
Bates  and  Peters  (1971)  in  the  development  of  rho-filtered  layergrams. 
They  formed  what  they  called  a  layergram  by  averaging  all  the  line  inte¬ 
gral  scans.  They  then  took  the  two-dimensional  Fourier  transform  of  this 
average,  weighted  the  transform  by  T  (p  in  their  notation),  and  inverted 
the  modified  transform  to  get  the  solution.  They  circumvented  the  slow¬ 
ness  of  a  two-dimensional  numerical  inversion  by  using  a  system  of  convex 
lenses  to  compute  the  transforms  optically. 


Convolution  method 

The  convolution  method  also  depends  on  the  convolution  theorem 
for  Fourier  transforms.  In  Gilbert’s  method  the  sum  of  the  line  integral 
scans  is  expressed  as  the  convolution  of  the  desired  solution  f  with  the 
inverse  transform  of  a  sampling  function  which  is  chosen  so  that  its 
inverse  transform  is  near  unity  in  the  domain  D.  In  the  convolution  or 
filtered  back  projection  method  the  scans  are  first  convolved  with  a 
filter  function  and  then  summed  to  yield  f.  The  basis  for  the  convolution 
technique  is  the  transformation  of  F(X,Y)  from  Cartesian  coordinates  X,Y 
to  polar  coordinates  T,0  which  results  in 


f(r) 


( T |  exp  (2iritT)  dT 


(25) 


where  t  *  x  cos9  +  y  sin0  and  j T |  is  the  Jacobian  of  the  transformation 
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The  equivalence  F(TC)  *  P(T,£)  allows  this  to  be  rewritten  as 


f(r> 


P(T,5)  |T|  exp(2ritT)  dT 


which  is,  by  virtue  of  the  convolution  theorem. 


f(r) 


p  *  h  d0 


p(t '  ,C)  h(t-t')  dt'. 


(26) 


(27) 


The  filter  function  h(t)  is  chosen  so  its  Fourier  transform  is 

H(T)  :  | T | . 


(28) 


Convolution  of  a  scan  with  such  a  filter  amounts  to  weighting  the  value 

of  the  integral  along  a  particular  line  L(t,0)  based  on  the  distance  of 

that  line  from  the  point  r  at  which  the  function  is  reconstructed.  The 

contributions  of  integrals  along  lines  near  r  are  emphasized  whereas 

these  at  a  distance  are  of  less  importance.  Generally  a  decrease  in 
-2 

the  amplitude  as  t  is  sufficient  to  provide  a  good  reconstruction. 

The  same  filter  function  is  convolved  with  the  scan  at  each  angle. 

Explicit  inversion  of  Eq.  (28)  to  obtain  h(t)  is  not  possible 
since  the  integral 


exp  (-2iritT)  dT 


(29) 


diverges;  therefore  an  approximation  of  the  integral  is  necessary.  One 

possible  approximation  for  the  integral  comes  from  limiting  the  range 

of  integration  to  a  finite  interval  [-Tq,  T  ].  The  resulting  integral 

will  converge  numerically  if  T  is  not  too  large.  Ramachandran  and 

o 
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Laks'nminarayanan  (197x/  evaluated  this  finite  integral  for  discrete  values 
of  t  separated  by  a  distance  b;  for  this  they  chose  To  =  (2b)  This 
discretization  reduces  the  inner  integral  in  Eq.  (27)  to  a  finite  sum. 

To  decrease  the  computation  time  further  they  assumed  h(t)  to  be  linear 
between  the  points  t=kb,  k=0,±l,  . ..,  ±b  .  This  assumption  allows  a 
linear  interpolation  between  the  values  t  *  kb  at  which  the  inner  sum 

in  Eq.  (27)  is  evaluated  and  the  values  t  =  x  cos9  +  y  sin0  needed  for 

the  reconstruction  of  f.  A  second  approximation  for  the  integral  (29) 
is  based  on  the  fact  that  F(T£),  and  hence  P(T,£)»  is  small  for  | T {  >  1. 

As  a  result  the  requirement  that  H(T)  ;  [ T J  may  be  weakened  for  | T (  >  1 

without  greatly  altering  the  inner  integral  in  Eq.  (26).  Tq  therefore  may 
be  taken  to  be  1  and  the  choice 

(  |T|  T  <  1 

H(T)  =  <  (30) 

(  0  T  5  1 

is  acceptable.  The  reconstruction  obtained  using  a  continuous  filter 

defined  in  this  way  (3racewell  and  Riddle  1967)  is  poor  because  h(t) 

decays  only  as  t  \  The  piecewise  linear  filter  of  Ramachandran  and 

-2 

Lakshminarayanan  which  decays  as  t  produces  better  reconstructicns . 

Shepp  and  Logan  (1974)  showed  that  the  assumption  of  pie-.ewise  linearity 
is  consistent  with  the  observation  that  H(T)  ;  |t|  for  |t|  <  1  and  there¬ 
fore  does  not  imply  a  loss  of  accuracy.  All  these  filters  have  a  large 
positive  amplitude  over  an  interval  centered  at  t=0  which  is  flanked  by 
negative  or  alternately  negative  and  positive  side  lobes  with  progressively 
smaller  amplitudes.  The  negative  side  lobes  tend  to  eliminate  the  star 
pattern  characteristic  of  the  back  projection  method,  as  illustrated  in 
Fig.  4.  In  this  diagram  the  solid  lines  mark  the  extent  of  the  positive 
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central  lobes  of  p^'  and  p^'.  The  heavy  shading  shows  regions  in  which 
this  overlapping  is  offset  by  negative  contributions  from  the  back  pro¬ 
jections  of  p^'  and  $2' •  This  reduction  in  the  intensity  of  the  overlap 
makes  it  easier  to  define  the  reconstructed  solution. 

a 

The  approximate  solution  f  obtained  using  the  convolution  method 
described  above  can  be  thought  of  as  a  blurred  version  of  the  true  solu¬ 
tion  f, 

f  -  f  *  *.  (31) 

Davison  and  Grunbaum  (1979)  chose  a  set  of  filter  functions  so 

that  the  reconstructed  function 

*  -  £  *  I  j,  Vi  •  <32> 

where  the  back  projection  operator  B ^  is  defined  by  (B^a  )(x,y)  =  a^(r*£  ), 
is  close  to  f  *  6  for  a  known  i{>.  The  point  response  function  <j>  is  chosen 
to  be  some  sort  of  identity  function  which  determines  the  resolution  of 
the  reconstruction.  This  approach  does  not  require  equally  spaced  angles 
as  do  the  convolution  algorithms  mentioned  previously.  Cho  et  at.  (1974) 
attempted  to  reduce  the  blurring  by  using  a  filter  based  on  a  geometrical 
correction  of  the  projection  data  in  the  spatial  domain. 


Series  expansion  methods 

In  the  series  expansion  approach  to  the  tomography  problem  the 
unknown  function  is  expressed  as  a  sum  of  known  orthogonal  basis  functions 
{bi(r)}> 


f  (r) 


(33) 


Applying  the  operator  L  results  in 


p(t.O  *  l  a  /  b  (r)  ds.  (34) 

i  1  A(t,e)  1  ' 

The  problem  is  then  reduced  to  finding  the  unknown  coefficients  a^.  If 
the  set  of  basis  functions  is  complete  the  reconstructed  solution  will 
be  unique.  Cormack  (1963)  expanded  the  solution  expressed  in  polar 
coordinates,  r,S  as  an  infinite  Fourier  series 


where 


f(r)  =  l  fn(r)  exp(in&) 

n=-c° 


/•2ir 

f n(t)  =  f(r)  exp(-inB)  dB  • 

0 


(35) 


(36) 


Under  the  assumption  that  the  f  's  are  bounded  and  piecewise  continuous, 

n 

a  similar  expansion  of  the  line  integral  scans  p(t,£)  results  in  a  unique 
relationship  between  the  Fourier  components  of  each  series.  Since  the 
set  of  Fourier  components  is  complete  the  reconstructed  solution  is 
unique.  In  a  later  paper,  Cormack  (1964)  extended  these  results  to 
include  the  entire  plane  and  considered  alternative  expansions  in  other 
orthogonal  functions.  For  a  set  of  discrete  angles  0^,  j=l,  ....  N, 

Eq.  (34)  is 

Pj(t)  =  l  a.  B±j (t)  j  =  1,  ...,  N  (37) 

where 

p,(t)  =  p(t,£  )  and  B.j(t)  =  I  b^r)  ds. 

■'Ut.e  ) 


Series  expansion  solutions  are  especially  useful  in  obtaining  con 
sistency  conditions  which  must  be  satisfied  if  a  set  of  functions  is  to 
represent  the  line  integral  scans  of  some  function.  Logan  and  Shepp 
(1975)  showed  that  the  minimum  norm  solution  to  the  discrete  problem 
when  the  angles  are  equally  spaced  is  of  the  form 

N 

f (r)  =  l  p . (r) 

j=l  J 


where  p ^ (r)  *  p(r*£j)  is  the  ridge  function  in  the  direction  0^  .  They 


sought  a  solution  by  expanding  p^  as  an  infinite  Fourier  sine  series. 


By  similarly  expanding  the  line  integral  scans  of  p^  they  obtained  cer¬ 


tain  consistency  conditions  which  must  be  satisfied  if  a  set  of  functions 


N 


{g(t,S^.)}j_^  is  to  correspond  to  a  unique  minimum  norm  solution.  These 


conditions  are  that:  1)  g(t,£.)  =  0  for  |t|  >1;  2)  the  integral 
i  J 

'22  -1/2 

[g(t,£  )]  (1-t  )  dt  is  finite;  and  3)  the  Fourier  components  of 


I 


g(t,i.)  satisfy  a  certain  relationship.  They  have  incorporated  the 


j 

symmetry  condition  g(t,£j)  =  g(-t,-C^)  by  restricting  the  range  of  angles 
to  [0,-).  Marr  (1974)  looked  for  a  least  squares  polynomial  solution 
to  the  fan  bean  tomography  problem.  The  consistency  conditions  imposed 
by  his  finite  Fourier  series  expansion  are  that  the  data  be  symmetric 
and  orthogonal  to  exp(in0).  Ludwig  (1966)  showed  that  the  consistency 
conditions  implied  by  an  expansion  in  spherical  harmonics  are  equivalent 
to  those  specified  above  for  the  general  solution  to  the  tomography 
problem.  In  a  preliminary  report,  Miller  (1978)  proposed  a  partial 
eigenfunction  expansion. 


Iterative  methods 


The  series  expansion  formulation  (37)  is  a  starting  point  for  the 
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iterative  and  matrix  inversion  reconstruction  techniques.  The  domain  D 
is  enclosed  by  a  qXq  grid  as  shown  in  Fig.  5  and  the  solution  f(r)  is 

2 

assumed  to  have  a  constant  mean  value  in  the  element  R^,  i*l,  ...,  q  . 
This  amounts  to  choosing  basis  functions 


bi(r) 


I1 

to 


i  -  1, 


2 

q  • 


(39) 


The  unknown  coefficients  a^  in  Eq.  (37)  are  then  the  f^'s.  The  projection 

data  are  regarded  as  sets  of  K  integrals  over  strips  S  ,  ,  kml,  ...,  K, 

j  k 

where  the  first  subscript  denotes  the  projection  angle  9  ^  ,  j=l,  ...»  N. 

The  total  number  of  data  is  then  KN  and  the  projections  may  be  renumbered 

KN 

to  form  the  set  {p^}^_^.  The  projection  operation  is  then 


i-1 


”ij£i 


i  - 1. 


,  KN 


(40) 


where  the  weight  function  w^_.  represents  the  fraction  of  the  iC^  element 

traversed  by  the  strip  S^.  Since  each  strip  intersects  only  a  few  of  the 

elements  the  matrix  W  composed  of  the  elements  w  is  quite  sparse. 

There  are  three  main  iterative  methods  which  have  been  applied  to 

£ 

the  tomography  problem.  They  differ  in  the  choice  of  a  correction  Af^ 
applied  to  f^,  the  iterate  of  the  solution,  to  obtain 


or 


fi+1  =  fi  +  ^fi 


f iZ+1  =  AfJ  •  fj 


i  -  1, 


(41) 


In  the  Iterative  Least  Squares  Technique  (ILST)  (Budinger  and  Cullberg 

£ 

1974b)  the  additive  correction  Af  is  chosen  so  as  to  minimize  the 

1  .  »  ; 
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squared  error  between  the  set  of  projections  of  the  solution  f  composed 
£  2 

of  all  elements  f^,  iml,  q  ,  and  the  observed  projections.  The 

Algebraic  Reconstruction  Technique  (ART)  (Gordon  et  dl.  1970)  uses  a 
projection  p^  to  correct  the  estimates  f^  for  elements  which  inter¬ 
sect  the  strip  S.,.  corresponding  to  p..  These  corrected  estimates  are 
j  k  j 

then  used  in  the  calculation  of  corrections  using  subsequent  projections. 
Hamaker  and  Solomon  (1978)  gave  a  rate  of  convergence  for  this  method 
which  depends  on  the  projection  angles.  The  Simultaneous  Iterative  Re¬ 
construction  Technique  (SIRT)  (Gilbert  1972a)  uses  all  the  projections 
to  correct  the  solution  in  each  element  R^.  ART  and  SIRT  have  both  addi¬ 
tive  and  multiplicative  formulations.  In  both  methods  the  additive  formu¬ 
lation  includes  a  positivity  constraint  on  the  f^'s.  The  multiplicative 
methods  will  produce  positive  f^'s  if  the  initial  guess  is  positive. 

Singular  Value  Decomposition 

Singular  value  decomposition  is  a  technique  for  finding  a  solution 

*  N 

f (;)  which  is  related  to  a  set  of  data  {y(Cj))j_^  by 

Y(C.)  =  f G^.Cj)  f(r)  dr  j  =  1,  ....  N  (42) 

•'v 

where  the  kernel  G(r,Cj)  describing  the  relevant  physical  process  in 
the  domain  V  is  non-symmetric.  In  the  notation  of  the  tomography  problem 
the  discrete  form  of  this  relation  is 


where  f  is  a  vector  composed  of  the  q  =  m  model  parameters  f^,  p  is  the 
data  vector  composed  of  the  KN  =  n  projections  p ^ ,  and  W  is  the  nXm  matrix 
whose  elements  are  the  weighting  factors  w^ .  If  W  were  non-singular. 
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symmetric  and  square  its  inverse  W  ^  would  exist  and  the  solution  to  Eq. 

(43)  would  be 

f  ■  W  ^  p.  (44) 

In  general  the  system  is  highly  underdetermined,  n  <<  m,  and  W  ^  does 
not  exist.  The  object  of  the  singular  value  decomposition  method  is 
the  construction  of  a  "pseudo-inverse"  for  a  rectangular  or ^ingular 
square  matrix.  The  following  discussion  is  based  on  that  of  Lanczos  (1961). 

The  eigenvalue  problem  for  the  symmetric  matrix  S  ~  1  reduces 

to 


*  S  '(?  o) 


W  v. 

'  -o 

T 

W  u . 


X  .  u . 
J  ~J 


X1  ?1 


j  *  1,  ...»  m 
i  -  1»  ...»  n 


(45) 


or  ecuivalently 


T  2 

W  W  u .  =  X  .  u . 

-  -  -l  l-i 


**  v  Yj  -  x32  Yj 


i  =  1,  ...»  n 

j  3  1»  •  •  • »  m 


(46) 


where  X^  =  X^  for  i= j .  The  values  X^  are  eigenvalues  of  the  symmetric 
T  T 

matrices  W  W  and  W  W;  the  X.'s  are  known  as  singular  values  of  W.  Let 
-  -  -  -  J 

SL  be  the  number  of  non-zero  eigenvalues.  Then  X^  ■  0  for  j  >  l  if  the 
singular  values  are  numbered  in  order  of  decreasing  magnitude.  Let  U 
and  V  be  unitary  matrices  whose  columns  are,  respectively,  the  eigen¬ 
vectors  u^  and  Vj  associated  with  the  l  non-zero  eigenvalues.  U  is  nXt 

and  V  is  mXS..  Let  U  and  V  denote,  respectively,  the  nX(n-L)  and 
-  -o  -o 

mX(m-Jl)  unitary  matrices  composed  of  the  remaining  eigenvectors.  The 
matrix  W  may  be  decomposed  into 


V  *  U  n  V 


(47) 
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where  Si  is  a  diagonal  matrix  whose  elements  are  the  singular  values  A^, 

J  ■  1»  1.  Although  the  operator  W  has  no  effect  on  elements  of 

Uq  and  Vq  these  subspaces  are  important  in  establishing  the  existence 

and  uniqueness  of  solutions  to  Wf  ■  g.  If  2  <  n,  a  solution  exists  if 

and  only  if  the  data  vector  p  is  orthogonal  to  all  solutions  of  the 

T 

homogeneous  equation  W  u^,  3  0,  that  is,  if  Uq  b  =  0.  If  2.  *  n,  then 
UQ  =  0  and  a  solution  will  always  exist.  If  n  <  £  <  m,  there  will  be 
m  -  2  solutions;  a  unique  solution  is  possible  only  if  2  =  m,  i.e.. 


V  =  0 . 
-o 


The  natural  inverse  or  the  generalized  inverse  (Penrose  1955) 


of  W  is 


+  -It 

w  «  v  n  u 


where  0.  ^  is  the  diagonal  matrix  composed  of  the  elements 

+  T 

j  =  1,  2.  The  product  W  W  =  U  U  may  be  thought  of  as  an  identity 

for  the  subspace  U  in  the  sense  that  multiplication  of  any  ueU  by  W  W+ 

reproduces  u.  If  W  W+  =  I  ,  then  U  =  0  and  the  consistency  condition 
-  *•  '  ~  n  ~o 

T 

tTo'  fc  =  0  is  identically  satisfied.  In  this  case  there  is  a  solution 
corresponding  to  any  data  vector  b.  The  product  W+  W  may  be  considered 


an  identity  for  V  in  a  similar  manner.  If  W  W  ■  I  ,  then  V  =  0  and  if 
*  *  %  -m  ~  o 

a  solution  exists  it  is  unique.  The  solution 

f  =  U+  p' 


is  seen  by  substitution  to  satisfy  Eq.  (43)  for  peU.  Thus  p  is  ortho¬ 
gonal  to  all  ueUq  and  satisfies  the  comparability  condition  for  the 
existence  of  a  solution.  Penrose  (1956)  shows  that  the  solution  (49) 
is  the  best  least  squares  solution  to  Eq.  (43)  in  the  sense  that 
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(fw  f  -  p||  >  ||w  f  -  p||  (50) 

or 

||  W  f  -  p | |  •  |[W  f  -  p | {  and  |  |  f  |  |  i  ||f||.  (51) 

where  ||x||  *  x  x,  for  any  other  solution  f.  . 

Budinger  and  Gullberg  (1974a)  used  the  generalized  inverse  to 
find  a  least  squares  solution  to  the  tonography  problem.  They  incorpo¬ 
rated  the  presence  of  errors  in  the  data  by  minimizing  instead  of 

T  T 

e  £,  where  £  *  Wf  -  p,  the  quantity  e  Pe  where  P  is  the  inverse  co- 

variance  matrix  for  the  data.  They  gave  the  solution 

f  -  (WT  P  W)"1  WT  P  p  (52) 

T 

or  if  U  P  W  is  singular,  as  is  usually  the  case, 

f  -  (P*  V)  +  f  p.  (53) 

Tewarson  (1972)  minimized  the  quadratic  form 

♦({)  *  £TP£  +  fTFf  (54) 

with  the  solution 

f  -  (WTPW  +  F)+  WTPp.  (55) 

The  matrices  ?  and  F  may  be  chosen  to  be  the  inverses  of  the  covariance 
matrices  of  the  data  and  the  model  parameters,  respectively,  or  they 
may  be  determined  using  the  regularization  method  of  Tikhonov  (Tikhonov 
and  Arsenin  1977). 

Evaluation  of  reconstruction  methods 


Reconstructions  based  on  projections  at  a  finite  number  of  angles 
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cannot  In  general  reproduce  exactly  the  source  function  and  approximate 
methods  for  reconstruction  are  generally  judged  by  the  closeness  of  the 
reconstructed  function  to  the  original  function.  One  technique  for 
evaluating  this  accuracy  is  to  measure  a  set  of  projections  of  a  physical 
object  with  known  properties,  called  a  phantom,  and  compare  the  object 
reconstructed  using  these  projections  with  the  original  object.  This 
evaluation  approach  is  deficient  in  that  the  inaccuracies  which  result 
may  be  due  not  to  the  reconstruction  method  but  to  experimental  and  com¬ 
putational  limitations.  Experimental  errors  are  eliminated  by  using 
instead  a  mathematical  phantom.  A  set  of  projections  of  a  known  function 
is  calculated  and  the  resulting  function  compared  with  the  source  function. 
Since  reconstructive  tomography  is  by  its  very  nature  a  computational 
phenomenon,  inaccuracies  due  to  round-off  and  other  purely  numerical 
errors  cannot  be  separated  from  inaccuracies  due  solely  to  the  recon¬ 
struction  method  itself.  The  accuracy  of  some  methods  may  be  character¬ 
ized  in  terms  of  an  error  bound  on  the  reconstructed  function  or,  in  the 
iterative  methods,  a  rate  of  convergence  of  the  algorithm.  Gordon  (1974) 
suggested  an  additional  assessment  of  the  accuracy  of  reconstructions 
using  iterative  techniques.  The  presence  of  some  feature  in  a  recon¬ 
structed  solution  may  be  tested  by  removing  the  feature  and  using  the 
altered  solution  as  an  initial  guess  for  a  second  application  of  the 
reconstruction  procedure.  If  the  feature  appears  again  in  the  second 
or  subsequent  reconstruction  it  is  not  likely  to  be  an  artifact  of  the 
reconstruction  method. 


CHAPTER  III 


STABILITY  OF  SOLUTIONS 


As  described  in  the  introduction  an  inverse  problem  is  well  posed 
if,  given  any  data  set,  a  solution  exists,  is  unique,  and  depends  continu¬ 
ously  on  the  data.  For  the  tomography  problem  Radon  (1917)  proved  that 
the  first  two  requirements  are  met  by  a  solution  f  subject  to  the  restric- 


be  finite,  but  he  did  not  examine  the  third  criterion.  Ludwig  (1966) 
showed  that  the  inversion  is  indeed  stable  within  the  space  of  rapidly 
decreasing  functions,  that  is-,  those  functions  f(r)  for  which  derivatives 
of  all  orders  exist,  are  continuous,  and  vanish  more  rapidly  than  r  , 
k  >  0,  as  r  +  *.  In  general,  however,  the  problem  of  inverting  the  line 
integral  data  to  obtain  f  is  not  well-posed. 

The  stability  of  the  discrete  problem 


W  f  =  p 


(56) 


with  a  finite  data  set  may  be  expressed  in  terms  of  the  singular  values 
of  the  matrix  W.  Since  the  problem  is  linear,  perturbations  6f  and  6p, 
not  necessarily  small,  irrespectively,  the  solution  and  the  data  are 
related  by 


W  <5f  -  6p 


(57) 
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T 

(Lanczoa  1961) .  The  decomposition  W  =  UOV  suggests  the  transformation 
6p  =  U<5p',  6f  *  VSf'  which  results  in 


fl  «'  -  6g'.  (58) 

Since  U  and  V  are  unitary  matrices  |5f'|  •  [  6  f  |  and  |§p'|  *  |  <5p  |  and  the 
variation  in  the  ith  component  of  the  solution  is 

5fi  =  Xi-1  6pi  *  (59) 


Thus  a  small  error  in  a  data  value  associated  with  a  small  singular 

value  can  cause  an  unacceptably  large  variation  in  the  solution. 

A  more  quantitative  measure  of  the  stability  of  a  discrete  inverse 

problem  is  the  variance  of  the  solution.  For  an  inverse  problem  to  be 

stable  the  variance  should  be  small.  The  variance  of  the  kth  component 

T 

of  the  solution  to  W  f  =  p  which  minimizes  f  f  is 


var  (fk) 


2 


a 


i 

l 


i-1 


(60) 


where  i  is  the  number  of  non-zero  singular  values  and  a  is  the  uncertainty 
in  the  data  due  to  measurement  errors  (Munk  and  Wunsch  1979) .  In  this 
formulation  the  errors  in  the  data  are  assumed  to  be  uncorrelated.  If 
the  errors  in  the  data  are  indeed  correlated  this  formulation  is  still 
valid  provided  certain  transformations  involving  the  covariance  matrices 
for  the  data  and  the  model  parameters  are  made  (Jackson  1972,  Munk  and 
Wunsch  1979).  If  one  or  more  of  the  singular  values  are  small  the  vari¬ 
ance  will  be  very  large.  One  way  to  reduce  the  variance  is  to  exclude 
very  small  singular  values  by  establishing  a  threshold  below  which  the 


singular  values  are  considered  zero.  The  ■'.atrix  V  would  then  be  composed 
of  the  eigenvectors  corresponding  to  singular  values  greater  than  the 
threshold  value.  The  threshold  may  be  chosen  so  the  uncertainty  in  the 
solution  is  sufficiently  small.  Unfortunately,  decreasing  the  number  of 
eigenvectors  in  V  also  decreases  the  resolution  in  the  solution.  A  unique, 
perfectly  resolved  solution  is  possible  only  when  none  of  the  singular 
values  vanish.  The  most  appropriate  compromise  in  this  trade-off  between 
minimizing  the  variance  and  maximizing  the  resolution  must  be  decided  for 
each  problem.  Davison  and  Griinbaum  (1979)  noted  a  similar  trade-off  between 
the  spatial  resolution  achieved  using  a  particular  convolution  algorithm 
and  the  closeness  of  the  point  response  function  for  that  algorithm  to 
a  specified  point  response  function. 

Smith  at  at,  (1977)  approached  the  problem  of  errors  in  the  data 
by  examining  whether  a  set  of  noisy  line  integrals  p(t,£j)  =  p^(t), 
j*i,  ...,  N,  can  actually  represent  the  line  integrals  of  some  function 
f.  They  proved  that  if  the  numbers 


-1 


(t)  t*  dt 


(61) 


are  the  values  at  the  points  5^  of  a  homogeneous  polynomial  of  degree 
k  £  N-2,  then  there  does  exist  a  function  f  whose  line  integral  scans 


in  the  directions  6^  are  the  given  data  p^(t).  This  conclusion  is  a  result 
of  the  form  of  the  operators  Ly  as  expressed  by  the  equivalence  (11). 


CHAPTER  IV 


PROSPECTS  FOR  FUTURE  RESEARCH 

Much  progress  has  been  made  on  the  tomography  problem,  but  some 
work  remains  to  be  done.  One  major  possibility  for  future  research  is 
the  incorporation  of  a  priori  information  into  the  problem.  There  are 
several  ways  in  which  this  may  be  done.  One  way  which  was  discussed 
previously  is  to  use  knowledge  of  the  symmetry  properties  of  the  function 
to  reduce  the  number  of  scans  needed  for  reconstruction.  DeRosier  and 
Kiug  (1968)  discussed  the  number  needed  in  electron  microscopy  to  re¬ 
construct  molecules  with  different  types  of  symmetry.  A  second  type  of 
z  priori  information  which  may  be  useful  is  a  bound  on  the  unknown 
function.  Miller  (1978)  proposed  using  the  boundedness  cf  the  norm 
of  f  as  a  stabilizing  constraint.  Davison  and  Grunbaum  (1979)  used 
this  bound  to  establish  an  error  bound  for  their  convolution  technique. 

A  third  possibility  is  the  incorporation  of  a  priori  knowledge  into  the 
construction  of  the  covariance  matrix  for  the  solution  (Munk  and  Wunsch 
1979) . 

Constraints  such  as  the  positivity  of  the  unknown  function  or 
the  total  mass  of  the  object  may  be  used  in  evaluating  the  validity  of 
a  reconstructed  solution.  A  solution  which  does  not  satisfy  these  con¬ 
straints  cannot  be  a  viable  solution.  So  far  the  positivity  constraint 
has  been  included  only  in  the  iterative  methods  (Gilbert  1972a,  Gordon 
et  at.  1970);  any  component  of  the  solution  which  turns  out  to  be 
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negative  is  set  to  zero.  The  problem  of  solving 

W  f  =  p  (62) 

subject  to  the  constraint 

f.  >  0  .  (63) 

l 

might  also  be  approached  using  techniques  in  linear  programming 
(Sabatier  1977a,  1977b)  in  which  the  constraint  is  incorporated  into 
the  formulation  of  the  problem. 
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-Convolution  or  filtered  back  projection  method.  Filtering  of  the  scans  reduces 
the  star  artifact  produced  by  the  back  projection  method.  See  text  for  explanation 


Spatial  discretization  of  reconstruction  domain  and  line  integral  data 
(after  Gordon  1974). 


